!pip install folium

import pandas as pd
import seaborn as sns
import folium
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.preprocessing import PolynomialFeatures
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import root_mean_squared_error, r2_score, mean_absolute_error, roc_curve, auc
from folium.plugins import HeatMap

# Setting up inline plotting
%matplotlib inline

data = pd.read_csv('dataset/housing_original.csv')
data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   longitude           20640 non-null  float64
 1   latitude            20640 non-null  float64
 2   housing_median_age  20640 non-null  float64
 3   total_rooms         20640 non-null  float64
 4   total_bedrooms      20433 non-null  float64
 5   population          20640 non-null  float64
 6   households          20640 non-null  float64
 7   median_income       20640 non-null  float64
 8   median_house_value  20640 non-null  float64
 9   ocean_proximity     20640 non-null  object 
dtypes: float64(9), object(1)
memory usage: 1.6+ MB

data.head()

data.describe()

# Drop rows with missing values
data.dropna(inplace=True)

# Drop duplicate rows
data.drop_duplicates(inplace=True)

# Define a function to remove outliers based on the IQR method
def remove_outliers(df, column):
    Q1 = df[column].quantile(0.25)
    Q3 = df[column].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    # Filter out rows with values outside the lower and upper bounds
    return df[(df[column] >= lower_bound) & (df[column] <= upper_bound)]

# Apply outlier removal for 'median_house_value'
data = remove_outliers(data, 'median_house_value').copy()

# Perform one-hot encoding for 'ocean_proximity'
data = pd.get_dummies(data, columns=['ocean_proximity'], drop_first=True)

# Check for any remaining missing values and confirm data structure
data.info()
data.head()
data.describe()

<class 'pandas.core.frame.DataFrame'>
Index: 19369 entries, 0 to 20639
Data columns (total 13 columns):
 #   Column                      Non-Null Count  Dtype  
---  ------                      --------------  -----  
 0   longitude                   19369 non-null  float64
 1   latitude                    19369 non-null  float64
 2   housing_median_age          19369 non-null  float64
 3   total_rooms                 19369 non-null  float64
 4   total_bedrooms              19369 non-null  float64
 5   population                  19369 non-null  float64
 6   households                  19369 non-null  float64
 7   median_income               19369 non-null  float64
 8   median_house_value          19369 non-null  float64
 9   ocean_proximity_INLAND      19369 non-null  bool   
 10  ocean_proximity_ISLAND      19369 non-null  bool   
 11  ocean_proximity_NEAR BAY    19369 non-null  bool   
 12  ocean_proximity_NEAR OCEAN  19369 non-null  bool   
dtypes: bool(4), float64(9)
memory usage: 1.6 MB

plt.style.use('ggplot')

plt.figure(figsize=(12, 8))
# Plotting the histogram of 'median_house_value' with more bins for granularity
sns.histplot(data['median_house_value'], kde=True, color='#3498db', bins=60)

# Adding vertical lines to show mean and median values for better insight
plt.axvline(x=data['median_house_value'].mean(), color='#e74c3c', linestyle='--', linewidth=2, label='Mean Value')
plt.axvline(x=data['median_house_value'].median(), color='#2ecc71', linestyle='-', linewidth=2, label='Median Value')
plt.title('Distribution of Median House Value in California', fontsize=18)
plt.xlabel('Median House Value (in USD)', fontsize=14)
plt.ylabel('Frequency', fontsize=14)
plt.legend(fontsize=12)
plt.grid(axis='y', linestyle='--', alpha=0.7)
mean_value = data['median_house_value'].mean()
median_value = data['median_house_value'].median()
plt.text(mean_value + 5000, plt.ylim()[1] * 0.7, f'Mean: ${mean_value:.0f}', color='red', fontsize=12, rotation=90)
plt.text(median_value + 5000, plt.ylim()[1] * 0.5, f'Median: ${median_value:.0f}', color='green', fontsize=12, rotation=90)
plt.show()
print("Figure 5.1.1: Distribution of Median House Value in California")

Figure 5.1.1: Distribution of Median House Value in California

# Distribution of Median Income Visualization
plt.figure(figsize=(12, 8))
sns.histplot(data['median_income'], kde=True, color='#3498db', bins=50)

# Adding vertical lines for the mean and median income for better insight
plt.axvline(x=data['median_income'].mean(), color='#e74c3c', linestyle='--', linewidth=2, label='Mean Income')
plt.axvline(x=data['median_income'].median(), color='#2ecc71', linestyle='-', linewidth=2, label='Median Income')
plt.title('Distribution of Median Income in California Housing Data', fontsize=18)
plt.xlabel('Median Income (in tens of thousands)', fontsize=14)
plt.ylabel('Frequency', fontsize=14)
plt.legend(fontsize=12)
plt.grid(axis='y', linestyle='--', alpha=0.7)

# Annotating the mean and median values
mean_income = data['median_income'].mean()
median_income = data['median_income'].median()
plt.text(mean_income + 0.1, plt.ylim()[1] * 0.7, f'Mean: ${mean_income:.2f}', color='red', fontsize=12, rotation=90)
plt.text(median_income + 0.1, plt.ylim()[1] * 0.5, f'Median: ${median_income:.2f}', color='blue', fontsize=12, rotation=90)
plt.show()
print("Figure 5.1.2: Distribution of Median Income in California Housing Data")

Figure 5.1.2: Distribution of Median Income in California Housing Data

plt.figure(figsize=(12, 8))
sns.scatterplot(x='median_income', y='median_house_value', data=data, alpha=0.6,
                hue='median_house_value', size='median_income', palette='viridis', sizes=(20, 200))

# Adding a regression line for better insight into the trend
sns.regplot(x='median_income', y='median_house_value', data=data, scatter=False, color='red', line_kws={'linestyle':'--'}, ci=None)
plt.title('Relationship between Median Income and Median House Value', fontsize=18)
plt.xlabel('Median Income (in tens of thousands)', fontsize=14)
plt.ylabel('Median House Value (in USD)', fontsize=14)
plt.legend(title='House Value', loc='upper left', fontsize=12, title_fontsize=14)
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
print("Figure 5.2.1: Relationship between Median Income and Median House Value")

Figure 5.2.1: Relationship between Median Income and Median House Value

# Create a figure with 1 row and 3 columns
fig, axs = plt.subplots(1, 3, figsize=(24, 8))

print("Figure 5.2.2: Relationship between Median Income and Median House Value by Ocean Proximity") 
# Plot 1: Median Income by Location
scatter1 = axs[0].scatter(data['longitude'], data['latitude'], c=data['median_income'], cmap='plasma', alpha=0.6, edgecolor='k', linewidth=0.5)
axs[0].set_xlabel('Longitude', fontsize=12)
axs[0].set_ylabel('Latitude', fontsize=12)
axs[0].set_title('Median Income by Location', fontsize=16)
fig.colorbar(scatter1, ax=axs[0], label='Income ($k)')

# Plot 2: Housing Median Age by Location
scatter2 = axs[1].scatter(data['longitude'], data['latitude'], c=data['housing_median_age'], cmap='cividis', alpha=0.8, edgecolor='k', linewidth=0.5)
axs[1].set_xlabel('Longitude', fontsize=12)
axs[1].set_ylabel('Latitude', fontsize=12)
axs[1].set_title('Housing Median Age by Location', fontsize=16)
fig.colorbar(scatter2, ax=axs[1], label='Housing Age (Years)')

# Plot 3: Median House Value by Location
scatter3 = axs[2].scatter(data['longitude'], data['latitude'], c=data['median_house_value'], cmap='magma', alpha=0.6, edgecolor='k', linewidth=0.5)
axs[2].set_xlabel('Longitude', fontsize=12)
axs[2].set_ylabel('Latitude', fontsize=12)
axs[2].set_title('Median House Value by Location', fontsize=16)
fig.colorbar(scatter3, ax=axs[2], label='House Value ($)')

plt.tight_layout(pad=3.0)
plt.subplots_adjust(top=0.92)
fig.suptitle('Geospatial Analysis of California Housing Data', fontsize=20, y=1.05)
plt.show()


# Creating a map centered around California
print("Figure 5.2.3: Heatmap of California Housing Data")
california_map = folium.Map(location=[34.986504,-118.716892], zoom_start = 6, min_zoom=4)
df_map = data[['latitude', 'longitude']]
heatmap_data = [[row['latitude'], row['longitude']] for index, row in df_map.iterrows()]
HeatMap(
    heatmap_data,
    radius=15,
    blur=10,
    max_zoom=1,
    min_opacity=0.3,
).add_to(california_map)
california_map

Figure 5.2.2: Relationship between Median Income and Median House Value by Ocean Proximity

Figure 5.2.3: Heatmap of California Housing Data

print("Figure 5.3.1: Correlation Heatmap of California Housing Dataset Features")
plt.figure(figsize=(14, 10))
sns.heatmap(
    data.corr(), 
    annot=True, 
    cmap='RdBu_r', # Reversed Red-Blue colormap to highlight both positive and negative correlations distinctly
    linewidths=1, 
    linecolor='white',
    annot_kws={'size': 10},
    fmt=".2f", 
    cbar_kws={'shrink': 0.8, 'aspect': 30, 'label': 'Correlation Coefficient'}
)
plt.title('Correlation Heatmap of California Housing Dataset Features', fontsize=18, pad=20)
plt.xticks(fontsize=12, rotation=45, ha='right')
plt.yticks(fontsize=12, rotation=0)
plt.tight_layout()
plt.show()

Figure 5.3.1: Correlation Heatmap of California Housing Dataset Features

data['ocean_proximity'] = data[['ocean_proximity_INLAND', 'ocean_proximity_ISLAND', 'ocean_proximity_NEAR BAY', 'ocean_proximity_NEAR OCEAN']].idxmax(axis=1).str.replace('ocean_proximity_', '')
print("Figure 5.4.1: Boxplot of Median House Value by Ocean Proximity")
sns.boxplot(x='ocean_proximity', y='median_house_value', data=data, hue='ocean_proximity', palette='Set3', width=0.6, fliersize=4, legend=False)
plt.title('Median House Value by Ocean Proximity', fontsize=20, pad=20)
plt.xlabel('Ocean Proximity', fontsize=15)
plt.ylabel('Median House Value (in USD)', fontsize=15)
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()

Figure 5.4.1: Boxplot of Median House Value by Ocean Proximity

plt.style.use('default')
selected_features = ['median_house_value', 'median_income', 'total_rooms', 'housing_median_age', 'population']
# Sampling the data for better readability and performance
sampled_data = data[selected_features].sample(5000, random_state=42)
print("Figure 5.5.1: Pair Plot of Selected Features from California Housing Data")
sns.pairplot(sampled_data, diag_kind='kde', kind='scatter', 
             plot_kws={'alpha': 0.5, 's': 15, 'edgecolor': 'k'}, diag_kws={'fill': True})
plt.suptitle('Pair Plot of Selected Features from California Housing Data', y=1.02, fontsize=18)
plt.show()

Figure 5.5.1: Pair Plot of Selected Features from California Housing Data

# Define the sample size and random state for splitting the data
sample_test_size = 0.8
sample_random_state = 42

# Selecting the feature and target for the first model: Predicting Median House Value using Median Income, Median Age, and Location
X = data[['longitude', 'latitude', 'housing_median_age','median_income']]
y = data['median_house_value'] #target variable

# Dictionary to store evaluation metrics
results = {
    'Model': [],
    'MSE': [],
    'R²': [],
    'MAE': []
}

# Helper function to evaluate and store results
def evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name):
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    mse = root_mean_squared_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    mae = mean_absolute_error(y_test, y_pred)
    if model_name in results['Model']:
        index = results['Model'].index(model_name)
        results['MSE'][index] = mse
        results['R²'][index] = r2
        results['MAE'][index] = mae
    else:
        results['Model'].append(model_name)
        results['MSE'].append(mse)
        results['R²'].append(r2)
        results['MAE'].append(mae)
    print(f"{str(model_no)}. {model_name} Linear Regression Model Evaluation:")
    print("Mean Squared Error (MSE):", mse)
    print("Mean Absolute Error (MAE):", mae)
    print("R-Squared (R²):", r2)
    return y_pred


# Helper function to plot residuals vs actual
def plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix):

    section = 'Figure 7.' + str(model_no)
    # Calculate residuals
    residuals = y_test - y_pred

    # Residual Distribution Plot    
    print(f"\n{section}.1: Residual Distribution for {model_name} Linear Regression ({title_suffix})")
    plt.style.use('ggplot')
    plt.figure(figsize=(14, 6))
    sns.histplot(residuals, kde=True, color='#3498db', bins=30)

    # Adding vertical lines for mean and median
    plt.axvline(residuals.mean(), color='#e74c3c', linestyle='--', linewidth=2, label='Mean')
    plt.axvline(residuals.median(), color='#2ecc71', linestyle='-', linewidth=2, label='Median')

    plt.title(f'Residual Distribution for {model_name} Linear Regression ({title_suffix})', fontsize=16)
    plt.xlabel('Residuals')
    plt.ylabel('Frequency')
    plt.legend()
    plt.grid(axis='y', linestyle='--', alpha=0.6)
    plt.show()

    # Predicted vs Actual Plot
    plt.figure(figsize=(16, 6))
    sc = plt.scatter(y_test, y_pred, alpha=0.6, c=y_pred, cmap='coolwarm', s=40)

    # Line of perfect prediction
    print(f"Figure {section}.2: Predicted vs. Actual Median House Value for {model_name} Linear Regression ({title_suffix})")
    plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', linewidth=2, label='Perfect Prediction')
    plt.xlabel('Actual Median House Value', fontsize=12)
    plt.ylabel('Predicted Median House Value', fontsize=12)
    plt.title(f'Predicted vs. Actual Median House Value ({model_name} Linear Regression)', fontsize=16)
    plt.colorbar(sc, label='Predicted House Value')
    plt.grid(True)
    plt.legend()
    plt.show()


# Heatmap for Predicted Median House Values
def plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix, poly_transformer=None):
    # Prepare data for heatmap
    section = 'Figure 7.' + str(model_no)

    print(f"{section}.3: Heatmap of {title_suffix} for {model_name} Linear Regression")

    # Merge y_test with X_test to get latitude and longitude
    y_test_df = pd.DataFrame(y_test).reset_index(drop=True)

    # Convert y_test to DataFrame and merge with X_test
    if poly_transformer:
        X_test_df = pd.DataFrame(X_test, columns=poly_transformer.get_feature_names_out(X.columns)).reset_index(drop=True)
        merged_df = pd.concat([X_test_df, y_test_df], axis=1)
        heat_data = [[row['latitude'], row['longitude'], prediction] for row, prediction in zip(merged_df.to_dict('records'), y_pred)]
    else:    
        y_test_df = pd.concat([y_test_df, X_test.reset_index(drop=True)], axis=1)
        heat_data = [[row['latitude'], row['longitude'], prediction] for row, prediction in zip(y_test_df.to_dict('records'), y_pred)]

    # Create base map
    map = folium.Map(location=[34.986504,-118.716892], zoom_start=6, min_zoom=4)

    # Add HeatMap layer
    HeatMap(heat_data, radius=8, blur=15, max_zoom=13).add_to(map)
    return map

model_no=1
model_name= 'Univariate'
title_suffix = 'Median Income vs. Median House Value'

# Splitting the dataset into training and testing sets
X_Uni = data[['median_income']]
X_train, X_test, y_train, y_test = train_test_split(X_Uni, y, test_size=sample_test_size, random_state=sample_random_state)

# Creating and training the Linear Regression model
model = LinearRegression()

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plat residuals vs. actual for the univariate model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

1. Univariate Linear Regression Model Evaluation:
Mean Squared Error (MSE): 72870.60208438961
Mean Absolute Error (MAE): 55706.32803939273
R-Squared (R²): 0.4153178441063634

Figure 7.1.1: Residual Distribution for Univariate Linear Regression (Median Income vs. Median House Value)

Figure Figure 7.1.2: Predicted vs. Actual Median House Value for Univariate Linear Regression (Median Income vs. Median House Value)

model_no=2
model_name= 'Bivariate'
title_suffix = 'Longitude & Latitude vs. Median House Value'

# Selecting features (longitude, latitude) and target (median house value)
X_loc = data[['longitude', 'latitude']]

# Splitting the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X_loc, y, test_size=sample_test_size, random_state=sample_random_state)

# Creating and training the Linear Regression model
model = LinearRegression()

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the bivariate model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

# Display the map
california_map = plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix)
california_map

2. Bivariate Linear Regression Model Evaluation:
Mean Squared Error (MSE): 82255.96212488494
Mean Absolute Error (MAE): 64215.74292257116
R-Squared (R²): 0.2550109369137644

Figure 7.2.1: Residual Distribution for Bivariate Linear Regression (Longitude & Latitude vs. Median House Value)

Figure Figure 7.2.2: Predicted vs. Actual Median House Value for Bivariate Linear Regression (Longitude & Latitude vs. Median House Value)

Figure 7.2.3: Heatmap of Longitude & Latitude vs. Median House Value for Bivariate Linear Regression

model_no = 3
model_name= 'Multivariate 1'
title_suffix = 'Income, Total Rooms, Population vs. Median House Value'

# define our predictor variables
X_multi = data[['median_income', 'total_rooms', 'population']]

# split data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X_multi, y, test_size=sample_test_size, random_state=sample_random_state)

# train the multiple linear regression model
model = LinearRegression()

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the Multivariate model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

3. Multivariate 1 Linear Regression Model Evaluation:
Mean Squared Error (MSE): 72922.80396725831
Mean Absolute Error (MAE): 55949.512409422096
R-Squared (R²): 0.41447985345442007

Figure 7.3.1: Residual Distribution for Multivariate 1 Linear Regression (Income, Total Rooms, Population vs. Median House Value)

Figure Figure 7.3.2: Predicted vs. Actual Median House Value for Multivariate 1 Linear Regression (Income, Total Rooms, Population vs. Median House Value)

model_no = 4
model_name= 'Multivariate 2'
title_suffix = 'Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value'


# training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=sample_test_size, random_state=sample_random_state)

# create a Linear Regression model
model = LinearRegression()

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the Multivariate model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

# Display the map
california_map = plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix)
california_map

4. Multivariate 2 Linear Regression Model Evaluation:
Mean Squared Error (MSE): 63197.01070376449
Mean Absolute Error (MAE): 47686.56049986861
R-Squared (R²): 0.560247622440619

Figure 7.4.1: Residual Distribution for Multivariate 2 Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure Figure 7.4.2: Predicted vs. Actual Median House Value for Multivariate 2 Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure 7.4.3: Heatmap of Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value for Multivariate 2 Linear Regression

model_no = 5
model_name = 'Ridge'
title_suffix = 'Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value'

# split data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=sample_test_size, random_state=sample_random_state)

# train the Ridge regression model
model = Ridge(alpha=1.0) 

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the Multivariate model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

# Display the map
california_map = plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix)
california_map

5. Ridge Linear Regression Model Evaluation:
Mean Squared Error (MSE): 63196.63009691232
Mean Absolute Error (MAE): 47686.2707734986
R-Squared (R²): 0.560252919281792

Figure 7.5.1: Residual Distribution for Ridge Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure Figure 7.5.2: Predicted vs. Actual Median House Value for Ridge Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure 7.5.3: Heatmap of Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value for Ridge Linear Regression

model_no = 6 
model_name = 'Lasso'
title_suffix = 'Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value'

# training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=sample_test_size, random_state=sample_random_state)

# Train the Lasso regression model
lasso_model = Lasso(alpha=1.0)

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the Lasso model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

# Display the map
california_map = plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix)
california_map

6. Lasso Linear Regression Model Evaluation:
Mean Squared Error (MSE): 63196.63009691232
Mean Absolute Error (MAE): 47686.2707734986
R-Squared (R²): 0.560252919281792

Figure 7.6.1: Residual Distribution for Lasso Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure Figure 7.6.2: Predicted vs. Actual Median House Value for Lasso Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure 7.6.3: Heatmap of Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value for Lasso Linear Regression

model_no = 7
model_name = 'Polynomial'
title_suffix = 'Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value'

# Transforming the feature set into polynomial features
poly_transformer = PolynomialFeatures(degree=2)  # Degree can be adjusted
X_poly_transformed = poly_transformer.fit_transform(X)

# Splitting the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X_poly_transformed, y, test_size=sample_test_size, random_state=sample_random_state)

# Creating and training the Polynomial Regression model
model = LinearRegression()

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the Lasso model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

# Display the map
california_map = plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix, poly_transformer)
california_map

7. Polynomial Linear Regression Model Evaluation:
Mean Squared Error (MSE): 60901.114569328405
Mean Absolute Error (MAE): 44511.59213809611
R-Squared (R²): 0.5916189271856771

Figure 7.7.1: Residual Distribution for Polynomial Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure Figure 7.7.2: Predicted vs. Actual Median House Value for Polynomial Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure 7.7.3: Heatmap of Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value for Polynomial Linear Regression

model_no = 8 
model_name = 'Random Forest'
title_suffix = 'Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value'

# Splitting the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=sample_test_size, random_state=sample_random_state)

# Creating and training the Random Forest Regressor model
model = RandomForestRegressor(n_estimators=100, random_state=sample_random_state)

# Evaluate the model
y_pred = evaluate_model(X_train, X_test, y_train, y_test, model, model_no, model_name)

# Plot residuals vs. actual for the Lasso model
plot_residuals_vs_actual(y_test, y_pred, model_no, model_name, title_suffix)

# Display the map
california_map = plot_geospatial_heatmap(X_test, y_test, y_pred, model_no, model_name, title_suffix)
california_map

8. Random Forest Linear Regression Model Evaluation:
Mean Squared Error (MSE): 47994.39438105085
Mean Absolute Error (MAE): 33100.54105962829
R-Squared (R²): 0.7463726081256825

Figure 7.8.1: Residual Distribution for Random Forest Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure Figure 7.8.2: Predicted vs. Actual Median House Value for Random Forest Linear Regression (Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value)

Figure 7.8.3: Heatmap of Median Income, Longitude, Latitude and Housing Median Age vs. Median House Value for Random Forest Linear Regression

# Convert results dictionary to DataFrame
results_df = pd.DataFrame(results)

# Displaying the results table for reference
print("\nFigure 8.1: Evaluation Metrics for models")
print(results_df)

print("\nFigure 8.2: Model Comparison for California Housing Data")

# Combined plot for Mean Squared Error (MSE), Mean Absolute Error (MAE), and R-Squared (R²)
plt.figure(figsize=(12, 8))
fig, ax1 = plt.subplots(figsize=(12, 8))

# Define color palette for a modern look
mse_color = '#66c2a5'   # Soft green for MSE
mae_color = '#fc8d62'   # Soft coral for MAE
r2_color = '#8da0cb'    # Soft blue for R²

# Plot for Mean Squared Error (MSE)
ax1.plot(results_df['Model'], results_df['MSE'], marker='o', linestyle='-', color=mse_color, linewidth=2, markersize=8, label='MSE')
ax1.set_xlabel("Model", fontsize=14)
ax1.set_ylabel("MSE", fontsize=14, color=mse_color)
ax1.tick_params(axis='y', labelcolor=mse_color)
for i, value in enumerate(results_df['MSE']):
    ax1.text(i, value, f"{value:.2e}", ha='center', va='bottom', fontsize=10, color=mse_color)
ax1.grid(True, linestyle='--', alpha=0.5)

# Plot for Mean Absolute Error (MAE)
ax2 = ax1.twinx()
ax2.plot(results_df['Model'], results_df['MAE'], marker='o', linestyle='-', color=mae_color, linewidth=2, markersize=8, label='MAE')
ax2.set_ylabel("MAE", fontsize=14, color=mae_color)
ax2.tick_params(axis='y', labelcolor=mae_color)
for i, value in enumerate(results_df['MAE']):
    ax2.text(i, value, f"{value:.2f}", ha='center', va='bottom', fontsize=10, color=mae_color)

# Plot for R-Squared (R²)
ax3 = ax1.twinx()
ax3.spines['right'].set_position(('outward', 60))
ax3.plot(results_df['Model'], results_df['R²'], marker='o', linestyle='-', color=r2_color, linewidth=2, markersize=8, label='R²')
ax3.set_ylabel("R²", fontsize=14, color=r2_color)
ax3.tick_params(axis='y', labelcolor=r2_color)
for i, value in enumerate(results_df['R²']):
    ax3.text(i, value, f"{value:.2f}", ha='center', va='bottom', fontsize=10, color=r2_color)

# Adding legends
fig.legend(loc='upper left', bbox_to_anchor=(0.1, 1), bbox_transform=ax1.transAxes, frameon=False)

# Title enhancement
plt.title("Model Comparison for California Housing Data", fontsize=16, weight='bold', color='#4a4a4a')

plt.tight_layout()
plt.show()

Figure 8.1: Evaluation Metrics for models
            Model           MSE        R²           MAE
0      Univariate  72870.602084  0.415318  55706.328039
1       Bivariate  82255.962125  0.255011  64215.742923
2  Multivariate 1  72922.803967  0.414480  55949.512409
3  Multivariate 2  63197.010704  0.560248  47686.560500
4           Ridge  63196.630097  0.560253  47686.270773
5           Lasso  63196.630097  0.560253  47686.270773
6      Polynomial  60901.114569  0.591619  44511.592138
7   Random Forest  47994.394381  0.746373  33100.541060

Figure 8.2: Model Comparison for California Housing Data

<Figure size 1200x800 with 0 Axes>

# Convert the regression problem to binary classification
# Threshold can be the median of median_house_value
median_value = data['median_house_value'].median()
y_binary = (data['median_house_value'] >= median_value).astype(int)  # 1 for high-value houses, 0 for low-value

# Prepare the features
X = data[['median_income', 'longitude', 'latitude', 'housing_median_age']]
X_train, X_test, y_train, y_test = train_test_split(X, y_binary, test_size=sample_test_size, random_state=sample_random_state)

# Define models
models = {
    "Univariate": LinearRegression(),
    "Bivariate": LinearRegression(),
    "Multivariate 1": LinearRegression(),
    "Multivariate 2": LinearRegression(),
    "Ridge": Ridge(alpha=1.0),
    "Lasso": Lasso(alpha=1.0),
    "Polynomial": PolynomialFeatures(degree=2),
    "Random Forest": RandomForestRegressor(n_estimators=100, random_state=42)
}

plt.figure(figsize=(16, 8))

# Loop through each model, fit, predict, and calculate ROC curve
for model_name, model in models.items():
    if model_name == "Polynomial":
        # Polynomial features transformation
        poly_transformer = PolynomialFeatures(degree=2)
        X_train_poly = poly_transformer.fit_transform(X_train)
        X_test_poly = poly_transformer.transform(X_test)
        poly_model = LinearRegression()
        poly_model.fit(X_train_poly, y_train)
        y_prob = poly_model.predict(X_test_poly)
    else:
        # Train other models
        model.fit(X_train, y_train)
        y_prob = model.predict(X_test)
    
    # Calculate ROC curve and AUC for the model
    fpr, tpr, _ = roc_curve(y_test, y_prob)
    roc_auc = auc(fpr, tpr)
    
    # Plot the ROC curve
    plt.plot(fpr, tpr, lw=2, label=f'{model_name} (AUC = {roc_auc:.2f})')

# Plot the diagonal line for random guessing
plt.plot([0, 1], [0, 1], 'k--', lw=2, label='Random Guessing (AUC = 0.50)')

# Adding plot details
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel("False Positive Rate", fontsize=12)
plt.ylabel("True Positive Rate", fontsize=12)
plt.title("Figure 8.2: Model Comparison for California Housing Data (ROC Curve)", fontsize=14)
plt.legend(loc="lower right", fontsize=10)
plt.grid(True, linestyle='--', alpha=0.7)
plt.show()

# Split the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=sample_test_size, random_state=sample_random_state)

# Function to perform hypothesis testing for OLS models
def hypothesis_test_ols(model):
    print(model.summary())  # This shows p-values, t-statistics, and more

# 1. Linear Regression Hypothesis Testing
X_train_const = sm.add_constant(X_train)
X_test_const = sm.add_constant(X_test)

linear_model = sm.OLS(y_train, X_train_const).fit()
hypothesis_test_ols(linear_model)  # Perform hypothesis testing for Linear Regression

# 2. Ridge Regression Hypothesis Testing (p-values not available directly in Ridge regression)
model = Ridge(alpha=1.0)
model.fit(X_train, y_train)

# Printing coefficients for Ridge, no p-values are available, so just the coefficients.
print("Ridge Regression Coefficients:", model.coef_)

# 3. Lasso Regression Hypothesis Testing (p-values not available directly in Lasso regression)
lasso_model = Lasso(alpha=1.0)
lasso_model.fit(X_train, y_train)

# Printing coefficients for Lasso
print("Lasso Regression Coefficients:", lasso_model.coef_)

# 4. Polynomial Regression Hypothesis Testing
poly_transformer = PolynomialFeatures(degree=2)
X_poly_train = poly_transformer.fit_transform(X_train)
X_poly_test = poly_transformer.transform(X_test)

poly_model = sm.OLS(y_train, X_poly_train).fit()
hypothesis_test_ols(poly_model)  # Perform hypothesis testing for Polynomial Regression

# 5. Random Forest Regression (Not suitable for hypothesis testing; feature importance can be used)
rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
rf_model.fit(X_train, y_train)

# Print feature importance
print("Random Forest Feature Importance:", rf_model.feature_importances_)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:     median_house_value   R-squared:                       0.558
Model:                            OLS   Adj. R-squared:                  0.558
Method:                 Least Squares   F-statistic:                     1221.
Date:                Fri, 18 Oct 2024   Prob (F-statistic):               0.00
Time:                        22:00:21   Log-Likelihood:                -48338.
No. Observations:                3873   AIC:                         9.669e+04
Df Residuals:                    3868   BIC:                         9.672e+04
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
======================================================================================
                         coef    std err          t      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------------
const              -3.623e+06   1.28e+05    -28.210      0.000   -3.87e+06   -3.37e+06
median_income       3.543e+04    696.151     50.888      0.000    3.41e+04    3.68e+04
longitude          -4.303e+04   1450.068    -29.677      0.000   -4.59e+04   -4.02e+04
latitude            -4.15e+04   1339.799    -30.972      0.000   -4.41e+04   -3.89e+04
housing_median_age   639.9831     87.514      7.313      0.000     468.405     811.561
==============================================================================
Omnibus:                      610.092   Durbin-Watson:                   1.993
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1163.007
Skew:                           0.976   Prob(JB):                    2.86e-253
Kurtosis:                       4.843   Cond. No.                     1.61e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.61e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Ridge Regression Coefficients: [ 35427.30136513 -42995.98199156 -41460.5514696     640.57925873]
Lasso Regression Coefficients: [ 35425.86090114 -43030.81672695 -41492.60712496    640.02879328]
                            OLS Regression Results                            
==============================================================================
Dep. Variable:     median_house_value   R-squared:                       0.594
Model:                            OLS   Adj. R-squared:                  0.592
Method:                 Least Squares   F-statistic:                     402.8
Date:                Fri, 18 Oct 2024   Prob (F-statistic):               0.00
Time:                        22:00:21   Log-Likelihood:                -48175.
No. Observations:                3873   AIC:                         9.638e+04
Df Residuals:                    3858   BIC:                         9.647e+04
Df Model:                          14                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       5.573e+07   1.04e+07      5.384      0.000    3.54e+07     7.6e+07
x1          -5.54e+05   9.25e+04     -5.991      0.000   -7.35e+05   -3.73e+05
x2          1.397e+06   2.25e+05      6.202      0.000    9.55e+05    1.84e+06
x3          1.533e+06   1.87e+05      8.216      0.000    1.17e+06     1.9e+06
x4         -6.103e+04   1.12e+04     -5.433      0.000   -8.31e+04    -3.9e+04
x5           123.5766    272.227      0.454      0.650    -410.147     657.300
x6         -6759.3123   1057.488     -6.392      0.000   -8832.601   -4686.024
x7         -6405.5299   1002.704     -6.388      0.000   -8371.411   -4439.649
x8           300.9656     58.779      5.120      0.000     185.725     416.206
x9          8744.9434   1248.565      7.004      0.000    6297.033    1.12e+04
x10         1.911e+04   2140.625      8.926      0.000    1.49e+04    2.33e+04
x11         -711.1416    128.805     -5.521      0.000    -963.675    -458.608
x12         1.034e+04   1021.434     10.120      0.000    8334.790    1.23e+04
x13         -724.9904    122.169     -5.934      0.000    -964.512    -485.468
x14           23.6474      6.332      3.735      0.000      11.234      36.061
==============================================================================
Omnibus:                      594.922   Durbin-Watson:                   2.000
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1301.433
Skew:                           0.902   Prob(JB):                    2.50e-283
Kurtosis:                       5.194   Cond. No.                     1.63e+08
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.63e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
Random Forest Feature Importance: [0.49234558 0.22531584 0.20984466 0.07249392]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity
0	-122.23	37.88	41.0	880.0	129.0	322.0	126.0	8.3252	452600.0	NEAR BAY
1	-122.22	37.86	21.0	7099.0	1106.0	2401.0	1138.0	8.3014	358500.0	NEAR BAY
2	-122.24	37.85	52.0	1467.0	190.0	496.0	177.0	7.2574	352100.0	NEAR BAY
3	-122.25	37.85	52.0	1274.0	235.0	558.0	219.0	5.6431	341300.0	NEAR BAY
4	-122.25	37.85	52.0	1627.0	280.0	565.0	259.0	3.8462	342200.0	NEAR BAY

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value
count	20640.000000	20640.000000	20640.000000	20640.000000	20433.000000	20640.000000	20640.000000	20640.000000	20640.000000
mean	-119.569704	35.631861	28.639486	2635.763081	537.870553	1425.476744	499.539680	3.870671	206855.816909
std	2.003532	2.135952	12.585558	2181.615252	421.385070	1132.462122	382.329753	1.899822	115395.615874
min	-124.350000	32.540000	1.000000	2.000000	1.000000	3.000000	1.000000	0.499900	14999.000000
25%	-121.800000	33.930000	18.000000	1447.750000	296.000000	787.000000	280.000000	2.563400	119600.000000
50%	-118.490000	34.260000	29.000000	2127.000000	435.000000	1166.000000	409.000000	3.534800	179700.000000
75%	-118.010000	37.710000	37.000000	3148.000000	647.000000	1725.000000	605.000000	4.743250	264725.000000
max	-114.310000	41.950000	52.000000	39320.000000	6445.000000	35682.000000	6082.000000	15.000100	500001.000000

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value
count	19369.000000	19369.000000	19369.000000	19369.000000	19369.000000	19369.000000	19369.000000	19369.000000	19369.000000
mean	-119.563902	35.655784	28.344158	2620.710930	539.893335	1442.285043	501.303991	3.665475	190802.064949
std	2.005895	2.151468	12.503931	2187.046669	422.650225	1145.780125	383.339200	1.556776	95404.934086
min	-124.350000	32.540000	1.000000	2.000000	2.000000	3.000000	2.000000	0.499900	14999.000000
25%	-121.760000	33.930000	18.000000	1440.000000	297.000000	798.000000	282.000000	2.522300	116100.000000
50%	-118.510000	34.270000	28.000000	2110.000000	437.000000	1181.000000	411.000000	3.442700	173200.000000
75%	-117.990000	37.730000	37.000000	3119.000000	648.000000	1746.000000	606.000000	4.572400	246400.000000
max	-114.310000	41.950000	52.000000	39320.000000	6445.000000	35682.000000	6082.000000	15.000100	482200.000000

California Housing Prices Statistical Analysis and Prediction¶

1. Introduction:¶

2. Objective:¶

3. Dataset:¶

Package Requirements:¶

Imports:¶

4. Data Cleaning:¶

5. Exploratory Data Analysis:¶

1. Distribution Analysis¶

2. Relationship Analysis¶

3. Correlation Analysis¶

4. Box Plot Analysis¶

5. Pair Plot Analysis¶

6. Exploratory Data Analysis Summary¶

7. Exploratory Data Analysis Conclusion¶

6. Model Selection:¶

Model Selection Criteria:¶

Planned Models for Selection:¶

Steps for Model Selection:¶

7. Model Analysis:¶

Model Analysis Criteria:¶

Sampling:¶

1. Univariate Linear Regression (Median House Value using Median Income)¶

2. Bivariate Linear Regression (Median House Value using Longitude and Latitude)¶

3. Multivariate Linear Regression (Median House Value using Median Income, Total Rooms, Population)¶

4. Multivariate Linear Regression (Median House Value using Median Income, Longitude, Latitude and Housing Median Age)¶

5. Ridge Regression (Median House Value using Median Income, Longitude, Latitude and Housing Median Age)¶

6. Lasso Regression (Median House Value using Median Income, Longitude, Latitude and Housing Median Age)¶

7. Polynomial Regression (Median House Value using Median Income, Longitude, Latitude and Housing Median Age)¶

8. Random Forest Regression (Median House Value using Median Income, Longitude, Latitude and Housing Median Age)¶

8. Comparative Analysis:¶

9. Hypothesis Validation and testing:¶

10. Conclusion:¶

11. References:¶